//----------------------------------*-C++-*----------------------------------//
/*!
 * \file   testSweepSource.cc
 * \author Jeremy Roberts
 * \date   Aug 31, 2011
 * \brief  Test of class SweepSource.
 * \note   Copyright (C) 2011 Jeremy Roberts. 
 */
//---------------------------------------------------------------------------//
// $Rev:: 161                                           $:Rev of last commit
// $Author:: j.alyn.roberts@gmail.com                   $:Author of last commit
// $Date:: 2011-10-24 11:45:46 +0000 (Mon, 24 Oct 2011) $:Date of last commit
//---------------------------------------------------------------------------//

#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "test/Testing.hh"
#include "test/TestInput.hh"

#include "angle/Definitions.hh"
#include "angle/MomentsIndex.hh"
#include "angle/Moment_to_Discrete.hh"

#include "material/BuiltInCrossSections.hh"
#include "material/Materials.hh"

#include "source/UniformSource.hh"

#include "tran/Moments_Field.hh"
#include "tran/CartMesh.hh"
#include "tran/Traits.hh"
#include "tran/Definitions.hh"
#include "tran/SweepSource.hh"

#include "../utilities/Constants.hh"
#include "../utilities/DBC.hh"
#include "../utilities/Soft_Equivalence.hh"
#include "../utilities/SP.hh"

// Equations
#include "tran/DD_1D_Equations.hh"
#include "tran/DD_2D_Equations.hh"
#include "tran/DD_3D_Equations.hh"

using namespace std;
using namespace slabtran;
using namespace util;
using namespace constants;

// typedefs
typedef SweepSource<DD_1D_Equations,_1D>::SP_sweepsource  SP_source_1d;
typedef SweepSource<DD_2D_Equations,_2D>::SP_sweepsource  SP_source_2d;
typedef SweepSource<DD_3D_Equations,_3D>::SP_sweepsource  SP_source_3d;

// 1d tests.
int test1d_construct(SP_source_1d &source);
int test1d_build_in_scatter_source(SP_source_1d &source);
int test1d_build_within_group_source(SP_source_1d &source);
int test1d_build_source(SP_source_1d &source);
// 2d tests
int test2d_construct(SP_source_2d &source);
int test2d_build_in_scatter_source(SP_source_2d &source);
int test2d_build_within_group_source(SP_source_2d &source);
int test2d_build_source(SP_source_2d &source);
// 3d tests
int test3d_construct(SP_source_3d &source);
int test3d_build_in_scatter_source(SP_source_3d &source);
int test3d_build_within_group_source(SP_source_3d &source);
int test3d_build_source(SP_source_3d &source);

int main(int argc, char *argv[])
{
  Require(argc==2);
  int test = atoi(argv[1]);
  cout << " ------ testSweepSource ------ " << endl;
  int test_value = 0;
  int number_of_tests = 13;

  // Declare sweep sources for tests.
  SP_source_1d s1; // 1d sweep source
  SP_source_2d s2; // 2d sweep source
  SP_source_3d s3; // 3d sweep source

  if ((test >=  1 && test <=  4) || test == number_of_tests) test_value += test1d_construct(s1);
  if (test ==  2 || test ==  4 || test == number_of_tests) test_value += test1d_build_in_scatter_source(s1);
  if (test ==  3 || test ==  4 || test == number_of_tests) test_value += test1d_build_within_group_source(s1);
  if (test ==  4 || test == number_of_tests              ) test_value += test1d_build_source(s1);

  if ((test >=  5 && test <=  8) || test == number_of_tests) test_value += test2d_construct(s2);
  if (test ==  6 || test ==  8 || test == number_of_tests) test_value += test2d_build_in_scatter_source(s2);
  if (test ==  7 || test ==  8 || test == number_of_tests) test_value += test2d_build_within_group_source(s2);
  if (test ==  8 || test == number_of_tests              ) test_value += test2d_build_source(s2);

  if ((test >=  9 && test <= 12) || test == number_of_tests) test_value += test3d_construct(s3);
  if (test == 10 || test == 12 || test == number_of_tests) test_value += test3d_build_in_scatter_source(s3);
  if (test == 11 || test == 12 || test == number_of_tests) test_value += test3d_build_within_group_source(s3);
  if (test == 12 || test == number_of_tests              ) test_value += test3d_build_source(s3);

  if (test > number_of_tests || test < 1)
    cout << " test " << test << " not available. " << endl;
  if (test_value)
    cout << " failed with value = " << test_value << endl;
  else
    cout << " passed." << endl;
  return test_value;
}

//===========================================================================//
// 1-D TESTS
//===========================================================================//
int test1d_construct(SP_source_1d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(6,          // number groups
                        1,          // legendre order
                        1,          // dimension
                        4,          // quad order
                        GAUSSLEG);  // quad type

  // Make a simple 1d mesh
  Mesh<_1D>::SP_mesh mesh;
  mesh = new Mesh<_1D>(10.0,        // width [cm]
                       8);          // meshes
  // Set all cells to material 0.
  Mesh<_1D>::Vec_Id matids(8, 0);
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_1D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_1D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_1D_Equations,_1D>::SP_state state;
  state = new State<DD_1D_Equations,_1D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_1D_Equations,_1D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS6G1M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Make a MomentsIndex
  MomentsIndex<_1D>::SP_momentsindex momentsindex;
  momentsindex = new MomentsIndex<_1D>(input->legendre_order());

  // Make an M operator
  Moment_to_Discrete<_1D>::SP_moment_to_discrete m_operator;
  m_operator = new Moment_to_Discrete<_1D>(input->legendre_order());
  m_operator->build(angularmesh);

  // Build the sweep source.
  source = new SweepSource<DD_1D_Equations,_1D>(state,
                                                mesh,
                                                angularmesh,
                                                materials,
                                                momentsindex,
                                                m_operator);
  // \note Here, we should add tests to ensure source is a valid
  //       sweep source.  Right now, there are no methods in
  //       \ref SweepSource that facilitate this, so we do nothing more.
  cout << " ... passed 1d SweepSource construct test. " << endl;
  return 0;
}

int test1d_build_in_scatter_source(SP_source_1d &source)
{
  Require(source);

  // Build the in scatter source.
  source->build_in_scatter_source(3);

  // Get the in scatter source.
  const SweepSource<DD_1D_Equations,_1D>::Moments_Field_t &s =
      source->in_scatter_source();
  TEST(s.legendre_order() == 1);

  // Set reference in-scatter source.
  double s_ref[] = {0.0, 0.0};
  // l = 0   phi_g0 * sig_{g->4,l=0}
  s_ref[0] = 1.0    * 0.42           +  // 2->4
             1.0    * 0.43           +  // 3->4
             1.0    * 0.45           ;  // 5->4
  // l = 1   phi_g1 * sig_{g->4,l=1}
  s_ref[1] = 1.0    * 0.042          +  // 2->4
             1.0    * 0.043          +  // 3->4
             1.0    * 0.045          ;  // 5->4

  // Check in-scatter source.
  SweepSource<DD_1D_Equations,_1D>::SP_mesh mesh =
      source->mesh();
  for (int i = 0; i < mesh->number_cells(); i++)
  {
    for (int l = 0; l <= s.legendre_order(); l++)
      TEST(soft_equiv(s(i, 0, l), s_ref[l])); // (cell, unknown, legendre)
  }

  cout << " ... passed 1d SweepSource build_in_scatter_source test. " << endl;
  return 0;
}

int test1d_build_within_group_source(SP_source_1d &source)
{
  Require(source);

  // Need the "working" flux for updating within group source, so
  //   get the state and use this group's flux already set to 1's
  State<DD_1D_Equations,_1D>::SP_state state = source->state();
  const State<DD_1D_Equations,_1D>::Moments_Field_t &phi = state->get_phi(3);

  // Build the within group scattering source.
  source->build_within_group_source(3, phi);

  // Get the within group scattering source.
  const SweepSource<DD_1D_Equations,_1D>::Moments_Field_t &s =
      source->within_group_source();
  TEST(s.legendre_order() == 1);

  // Set reference within-group scatter source.
  double s_ref[] = {0.0, 0.0};
  // l = 0   phi_g0 * sig_{g->4,l=0}
  s_ref[0] = 1.0    * 0.44           ;  // 4->4
  // l = 1   phi_g1 * sig_{g->4,l=1}
  s_ref[1] = 1.0    * 0.044          ;  // 4->4

  // Check within-group scatter source.
  SweepSource<DD_1D_Equations,_1D>::SP_mesh mesh =
      source->mesh();
  for (int i = 0; i < mesh->number_cells(); i++)
  {
    for (int l = 0; l <= s.legendre_order(); l++)
      TEST(soft_equiv(s(i, 0, l), s_ref[l])); // (cell, unknown, legendre)
  }

  cout << " ... passed 1d SweepSource build_within_group_source test. " << endl;
  return 0;
}

int test1d_build_source(SP_source_1d &source)
{
  Require(source);

  // Build the sweep source for octant 0 and angle 0.
  source->build_source(3, 0, 0);

  // Get the sweep source.
  const SweepSource<DD_1D_Equations,_1D>::SweepSourceField &s =
      source->source();
  TEST(s.size() == (source->mesh())->number_cells());

  // Set reference source.
  //   First angle in the quadrature.
  double mu = 0.8611363115940525752239465;
  //   In-scatter component.
  //             1/2 * phi_g0 * sig_{g->4,l=0} * P_0(mu)    3/2 * phi_g1 * sig_{g->4,l=1} * P_1(mu)
  double s_ref = 0.5 * 1.0    * 0.42           * 1.0      + 1.5 * 1.0    * 0.042          * mu      +  // 2->4
                 0.5 * 1.0    * 0.43           * 1.0      + 1.5 * 1.0    * 0.043          * mu      +  // 3->4
                 0.5 * 1.0    * 0.45           * 1.0      + 1.5 * 1.0    * 0.045          * mu      ;  // 5->4
  //   Within-group component.
  s_ref       += 0.5 * 1.0    * 0.44           * 1.0      + 1.5 * 1.0    * 0.044          * mu      ;  // 4->4

  // Check sweep source.
  for (int i = 0; i < s.size(); i++)
    TEST(soft_equiv(s[i], s_ref));

  // Create an external source and add it to source.
  UniformSource<_1D>::SP_source qe;
  qe = new UniformSource<_1D>(6, 1.0);
  source->set_external_source(qe);

  // Update the reference source.  (Remember it's divided by two).
  s_ref       += 0.5;

  // Rebuild sweep source and check.
  source->build_source(3, 0, 0);
  for (int i = 0; i < s.size(); i++)
    TEST(soft_equiv(s[i], s_ref));

  cout << " ... passed 1d SweepSource build_source test. " << endl;
  return 0;
}

//===========================================================================//
// 2-D TESTS
//===========================================================================//

int test2d_construct(SP_source_2d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(6,          // number groups
                        1,          // legendre order
                        2,          // dimension
                        4,          // quad order
                        LEVELSYM);  // quad type

  // Make a simple 2d mesh
  Mesh<_2D>::SP_mesh mesh;
  mesh = new Mesh<_2D>(10.0,        // width [cm]
                       8);          // meshes
  // Set all cells to material 0.
  Mesh<_2D>::Vec_Id matids(mesh->number_cells(), 0);
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_2D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_2D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_2D_Equations,_2D>::SP_state state;
  state = new State<DD_2D_Equations,_2D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_2D_Equations,_2D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS6G1M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Make a MomentsIndex
  MomentsIndex<_2D>::SP_momentsindex momentsindex;
  momentsindex = new MomentsIndex<_2D>(input->legendre_order());

  // Make an M operator
  Moment_to_Discrete<_2D>::SP_moment_to_discrete m_operator;
  m_operator = new Moment_to_Discrete<_2D>(input->legendre_order());
  m_operator->build(angularmesh);

  // Build the sweep source.
  source = new SweepSource<DD_2D_Equations,_2D>(state,
                                                mesh,
                                                angularmesh,
                                                materials,
                                                momentsindex,
                                                m_operator);
  // \note Here, we should add tests to ensure source is a valid
  //       sweep source.  Right now, there are no methods in
  //       \ref SweepSource that facilitate this, so we do nothing more.
  cout << " ... passed 2d SweepSource construct test. " << endl;
  return 0;
}

int test2d_build_in_scatter_source(SP_source_2d &source)
{
  Require(source);

  // Build the in scatter source.
  source->build_in_scatter_source(3);

  // Get the in scatter source.
  const SweepSource<DD_2D_Equations,_2D>::Moments_Field_t &s =
      source->in_scatter_source();
  TEST(s.legendre_order() == 1);

  // Set reference in-scatter source.
  double s_ref[] = {0.0, 0.0, 0.0};
  // l,m = 0,0   phi^0_0g   * sig_{g->4,l=0}
  s_ref[0]     = 1.0        * 0.42           +  // 2->4
                 1.0        * 0.43           +  // 3->4
                 1.0        * 0.45           ;  // 5->4
  // l,m = 1,-1  phi^-1_1g * sig_{g->4,l=1}
  s_ref[1]     = 1.0        * 0.042          +  // 2->4
                 1.0        * 0.043          +  // 3->4
                 1.0        * 0.045          ;  // 5->4
  // l,m = 1,1   phi^+1_1g * sig_{g->4,l=1}
  s_ref[2]     = 1.0        * 0.042          +  // 2->4
                 1.0        * 0.043          +  // 3->4
                 1.0        * 0.045          ;  // 5->4
  // Check in-scatter source.
  SweepSource<DD_2D_Equations,_2D>::SP_mesh mesh =
      source->mesh();
  SweepSource<DD_2D_Equations,_2D>::SP_momentsindex momentsindex =
      source->momentsindex();
  for (int i = 0; i < mesh->number_cells(); i++) // cardinal cell index
  {
    int j = 0; // counter for reference source
    for (int l = 0; l <= s.legendre_order(); l++)
    {
      // Get m index.
      const vector<int> &m = momentsindex->m_index(l);
      for (int i = 0; i < m.size(); i++)
        TEST(soft_equiv(s(i, 0, l, m[i]), s_ref[j++])); // (cell, unknown, l, m)
    }
  }

  cout << " ... passed 2d SweepSource build_in_scatter_source test. " << endl;
  return 0;
}

int test2d_build_within_group_source(SP_source_2d &source)
{
  Require(source);

  // Need the "working" flux for updating within group source, so
  // get the state and use this group's flux already set to 1's
  State<DD_2D_Equations,_2D>::SP_state state = source->state();
  const State<DD_2D_Equations,_2D>::Moments_Field_t &phi = state->get_phi(3);

  // Build the within group scattering source.
  source->build_within_group_source(3, phi);

  // Get the within group scattering source.
  const SweepSource<DD_2D_Equations,_2D>::Moments_Field_t &s =
      source->within_group_source();
  TEST(s.legendre_order() == 1);

  // Set reference within-group scatter source.
  double s_ref[] = {0.0, 0.0, 0.0};
  // l,m = 0,0   phi^0_0g   * sig_{g->4,l=0}
  s_ref[0]     = 1.0        * 0.44           ;  // 4->4
  // l,m = 1,-1  phi^-1_1g * sig_{g->4,l=1}
  s_ref[1]     = 1.0        * 0.044          ;  // 4->4
  // l,m = 1,1   phi^-1_1g * sig_{g->4,l=1}
  s_ref[2]     = 1.0        * 0.044          ;  // 4->4

  // Check within-group scatter source.
  SweepSource<DD_2D_Equations,_2D>::SP_mesh mesh =
      source->mesh();
  SweepSource<DD_2D_Equations,_2D>::SP_momentsindex momentsindex =
      source->momentsindex();
  for (int i = 0; i < mesh->number_cells(); i++) // cardinal cell index
  {
    int j = 0; // counter for reference source
    for (int l = 0; l <= s.legendre_order(); l++)
    {
      // Get m index.
      const vector<int> &m = momentsindex->m_index(l);
      for (int i = 0; i < m.size(); i++)
        TEST(soft_equiv(s(i, 0, l, m[i]), s_ref[j++])); // (cell, unknown, l, m)
    }
  }

  cout << " ... passed 2d SweepSource build_within_group_source test. " << endl;
  return 0;
}

int test2d_build_source(SP_source_2d &source)
{
  Require(source);

  // Build the sweep source for octant 0 and angle 0.
  source->build_source(3, 0, 0);

  // Get the sweep source.
  const SweepSource<DD_2D_Equations,_2D>::SweepSourceField &s =
      source->source();
  TEST(s.size() == (source->mesh())->number_cells());

  // Set reference source.
  //   First angle in the quadrature.
  double mu     = 0.3500211745815407;
  double eta    = 0.3500211745815407;
  //   In-scatter component.  Divide by 4*pi after.
  //             1.0 * p[0, 0] * s_{g->4,l=0} * Y[0, 0]  + 3.0 * p[1,-1] * s_{g->4,l=1} * Y[1,-1]  + 3.0 * p[1, 1] * s_{g->4,l=1} * Y[1, 1]
  double s_ref = 1.0 * 1.0    * 0.42          * 1.0      + 3.0 * 1.0     * 0.042        * mu       + 3.0 * 1.0     * 0.042        * eta   + // 2->4
                 1.0 * 1.0    * 0.43          * 1.0      + 3.0 * 1.0     * 0.043        * mu       + 3.0 * 1.0     * 0.043        * eta   + // 3->4
                 1.0 * 1.0    * 0.45          * 1.0      + 3.0 * 1.0     * 0.045        * mu       + 3.0 * 1.0     * 0.045        * eta   ; // 5->4
  //   Within-group component.
  s_ref       += 1.0 * 1.0    * 0.44          * 1.0      + 3.0 * 1.0     * 0.044        * mu       + 3.0 * 1.0     * 0.044        * eta   ; // 4->4
  //   Normalize by 4*pi.
  s_ref       *= inv_four_pi;

  // Check sweep source.
  for (int i = 0; i < s.size(); i++)
    TEST(soft_equiv(s[i], s_ref));

  // Create an external source and add it to source.
  UniformSource<_2D>::SP_source qe;
  qe = new UniformSource<_2D>(6, 1.0);
  source->set_external_source(qe);

  // Update the reference source.  (Remember it's divided by 4*pi).
  s_ref       += inv_four_pi;

  // Rebuild sweep source and check.
  source->build_source(3, 0, 0);
  for (int i = 0; i < s.size(); i++)
    TEST(soft_equiv(s[i], s_ref));

  cout << " ... passed 2d SweepSource build_source test. " << endl;
  return 0;
}


//===========================================================================//
// 3-D TESTS
//===========================================================================//

int test3d_construct(SP_source_3d &source)
{
  // Make test input.
  TestInput::SP_input input;
  input = new TestInput(6,          // number groups
                        1,          // legendre order
                        3,          // dimension
                        4,          // quad order
                        LEVELSYM);  // quad type

  // Make a simple 3d mesh
  Mesh<_3D>::SP_mesh mesh;
  mesh = new Mesh<_3D>(10.0,        // width [cm]
                       8);          // meshes
  // Set all cells to material 0.
  Mesh<_3D>::Vec_Id matids(mesh->number_cells(), 0);
  mesh->set_cell_id(MATERIAL_ID, matids);

  // Make a simple angular mesh.
  AngularMesh<_3D>::SP_angularmesh angularmesh;
  angularmesh = new AngularMesh<_3D>(input->quadrature_order(),
                                     input->quadrature_type());

  // Make a state and set the moments to unity.
  State<DD_3D_Equations,_3D>::SP_state state;
  state = new State<DD_3D_Equations,_3D>(input, mesh, angularmesh);
  for (int g = 0; g < state->number_groups(); g++)
  {
    State<DD_3D_Equations,_3D>::Moments_Field_t &phi = state->get_phi(g);
    for (int i = 0; i < phi.size(); i++)
      phi[i] = 1.0;
  }

  // Use a built-in library.
  BuiltInCrossSections xs(BuiltInCrossSections::XS6G1M);
  Materials::SP_materials materials;
  materials = xs.build();

  // Make a MomentsIndex
  MomentsIndex<_3D>::SP_momentsindex momentsindex;
  momentsindex = new MomentsIndex<_3D>(input->legendre_order());

  // Make an M operator
  Moment_to_Discrete<_3D>::SP_moment_to_discrete m_operator;
  m_operator = new Moment_to_Discrete<_3D>(input->legendre_order());
  m_operator->build(angularmesh);

  // Build the sweep source.
  source = new SweepSource<DD_3D_Equations,_3D>(state,
                                                mesh,
                                                angularmesh,
                                                materials,
                                                momentsindex,
                                                m_operator);

  // \note Here, we should add tests to ensure source is a valid
  //       sweep source.  Right now, there are no methods in
  //       \ref SweepSource that facilitate this, so we do nothing more.
  cout << " ... passed 3d SweepSource construct test. " << endl;
  return 0;
}

int test3d_build_in_scatter_source(SP_source_3d &source)
{
  Require(source);

  // Build the in scatter source.
  source->build_in_scatter_source(3);

  // Get the in scatter source.
  const SweepSource<DD_3D_Equations,_3D>::Moments_Field_t &s =
      source->in_scatter_source();
  TEST(s.legendre_order() == 1);

  // Set reference in-scatter source.
  double s_ref[] = {0.0, 0.0, 0.0, 0.0};
  // l,m = 0,0   phi^0_0g   * sig_{g->4,l=0}
  s_ref[0]     = 1.0        * 0.42           +  // 2->4
                 1.0        * 0.43           +  // 3->4
                 1.0        * 0.45           ;  // 5->4
  // l,m = 1,-1  phi^-1_1g  * sig_{g->4,l=1}
  s_ref[1]     = 1.0        * 0.042          +  // 2->4
                 1.0        * 0.043          +  // 3->4
                 1.0        * 0.045          ;  // 5->4
  // l,m = 1,0   phi^0_1g   * sig_{g->4,l=1}
  s_ref[2]     = 1.0        * 0.042          +  // 2->4
                 1.0        * 0.043          +  // 3->4
                 1.0        * 0.045          ;  // 5->4
  // l,m = 1,1   phi^-1_1g  * sig_{g->4,l=1}
  s_ref[3]     = 1.0        * 0.042          +  // 2->4
                 1.0        * 0.043          +  // 3->4
                 1.0        * 0.045          ;  // 5->4

  // Check in-scatter source.
  SweepSource<DD_3D_Equations,_3D>::SP_mesh mesh =
      source->mesh();
  SweepSource<DD_3D_Equations,_3D>::SP_momentsindex momentsindex =
      source->momentsindex();
  TEST(s.size() == 4*mesh->number_cells());
  for (int i = 0; i < mesh->number_cells(); i++) // cardinal cell index
  {
    int j = 0; // counter for reference source
    for (int l = 0; l <= s.legendre_order(); l++)
    {
      // Get m index.
      const vector<int> &m = momentsindex->m_index(l);
      for (int i = 0; i < m.size(); i++)
        TEST(soft_equiv(s(i, 0, l, m[i]), s_ref[j++])); // (cell, unknown, l, m)
    }
  }

  cout << " ... passed 3d SweepSource build_in_scatter_source test. " << endl;
  return 0;
}

int test3d_build_within_group_source(SP_source_3d &source)
{
  Require(source);
  // Need the "working" flux for updating within group source, so
  // get the state and use this group's flux already set to 1's
  State<DD_3D_Equations,_3D>::SP_state state = source->state();
  const State<DD_3D_Equations,_3D>::Moments_Field_t &phi = state->get_phi(3);

  // Build the within group scattering source.
  source->build_within_group_source(3, phi);

  // Get the within group scattering source.
  const SweepSource<DD_3D_Equations,_3D>::Moments_Field_t &s =
      source->within_group_source();
  TEST(s.legendre_order() == 1);

  // Set reference within-group scatter source.
  double s_ref[] = {0.0, 0.0, 0.0};
  // l,m = 0,0   phi^0_0g   * sig_{g->4,l=0}
  s_ref[0]     = 1.0        * 0.44           ;  // 4->4
  // l,m = 1,-1  phi^-1_1g  * sig_{g->4,l=1}
  s_ref[1]     = 1.0        * 0.044          ;  // 4->4
  // l,m = 1,0   phi^0_1g   * sig_{g->4,l=1}
  s_ref[2]     = 1.0        * 0.044          ;  // 4->4
  // l,m = 1,1   phi^-1_1g  * sig_{g->4,l=1}
  s_ref[3]     = 1.0        * 0.044          ;  // 4->4

  // Check within-group scatter source.
  SweepSource<DD_3D_Equations,_3D>::SP_mesh mesh =
      source->mesh();
  SweepSource<DD_3D_Equations,_3D>::SP_momentsindex momentsindex =
      source->momentsindex();
  for (int i = 0; i < mesh->number_cells(); i++) // cardinal cell index
  {
    int j = 0; // counter for reference source
    for (int l = 0; l <= s.legendre_order(); l++)
    {
      // Get m index.
      const vector<int> &m = momentsindex->m_index(l);
      for (int i = 0; i < m.size(); i++)
        TEST(soft_equiv(s(i, 0, l, m[i]), s_ref[j++])); // (cell, unknown, l, m)
    }
  }

  cout << " ... passed 3d SweepSource build_within_group_source test. " << endl;
  return 0;
}

int test3d_build_source(SP_source_3d &source)
{
  Require(source);

  // Build the sweep source for octant 0 and angle 0.
  source->build_source(3, 0, 0);

  // Get the sweep source.
  const SweepSource<DD_3D_Equations,_3D>::SweepSourceField &s =
      source->source();
  TEST(s.size() == (source->mesh())->number_cells());

  // Set reference source.
  //   First angle in the quadrature.
  double mu     = 0.3500211745815407;
  double eta    = 0.3500211745815407;
  double xi     = 0.8688903007222012;
  //   In-scatter component.  Divide by 4*pi after.
  //             1.0 * p[0, 0] * s_{g->4,l=0} * Y[0, 0]  + 3.0 * p[1,-1] * s_{g->4,l=1} * Y[1,-1]  + 3.0 * p[1, 0] * s_{g->4,l=1} * Y[1, 0]  + 3.0 * p[1, 1] * s_{g->4,l=1} * Y[1, 1]
  double s_ref = 1.0 * 1.0    * 0.42          * 1.0      + 3.0 * 1.0     * 0.042        * mu       + 3.0 * 1.0     * 0.042        * xi       + 3.0 * 1.0     * 0.042        * eta   + // 2->4
                 1.0 * 1.0    * 0.43          * 1.0      + 3.0 * 1.0     * 0.043        * mu       + 3.0 * 1.0     * 0.043        * xi       + 3.0 * 1.0     * 0.043        * eta   + // 3->4
                 1.0 * 1.0    * 0.45          * 1.0      + 3.0 * 1.0     * 0.045        * mu       + 3.0 * 1.0     * 0.045        * xi       + 3.0 * 1.0     * 0.045        * eta   ; // 5->4
  //   Within-group component.
  s_ref       += 1.0 * 1.0    * 0.44          * 1.0      + 3.0 * 1.0     * 0.044        * mu       + 3.0 * 1.0     * 0.044        * xi       + 3.0 * 1.0     * 0.044        * eta   ; // 4->4
  //   Normalize by 4*pi.
  s_ref       *= inv_four_pi;

  // Check sweep source.
  for (int i = 0; i < s.size(); i++)
  {
    TEST(soft_equiv(s[i], s_ref));
  }

  // Create an external source and add it to source.
  UniformSource<_3D>::SP_source qe;
  qe = new UniformSource<_3D>(6, 1.0);
  source->set_external_source(qe);

  // Update the reference source.  (Remember it's divided by 4*pi).
  s_ref       += inv_four_pi;

  // Rebuild sweep source and check.
  source->build_source(3, 0, 0);
  for (int i = 0; i < s.size(); i++)
    TEST(soft_equiv(s[i], s_ref));

  cout << " ... passed 3d SweepSource build_source test. " << endl;
  return 0;
}

//---------------------------------------------------------------------------//
//              end of testSweepSource.cc
//---------------------------------------------------------------------------//
